{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d83bafb9",
   "metadata": {},
   "source": [
    "# 习题\n",
    "## 习题10.1\n",
    "![image.png](./images/exercise1.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "318a5073",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "初始化后向概率：\n",
      "beta4(1) = 1\n",
      "beta4(2) = 1\n",
      "beta4(3) = 1\n",
      "\n",
      "从时刻T-1到1观测序列的后向概率：\n",
      "beta3(1) = sum[a1j * bj(o4) * beta4(j)] = (0.50 * 0.50 * 1.00 + 0.20 * 0.60 * 1.00 + 0.30 * 0.30 * 1.00 + 0) = 0.460\n",
      "beta3(2) = sum[a2j * bj(o4) * beta4(j)] = (0.30 * 0.50 * 1.00 + 0.50 * 0.60 * 1.00 + 0.20 * 0.30 * 1.00 + 0) = 0.510\n",
      "beta3(3) = sum[a3j * bj(o4) * beta4(j)] = (0.20 * 0.50 * 1.00 + 0.30 * 0.60 * 1.00 + 0.50 * 0.30 * 1.00 + 0) = 0.430\n",
      "beta2(1) = sum[a1j * bj(o3) * beta3(j)] = (0.50 * 0.50 * 0.46 + 0.20 * 0.40 * 0.51 + 0.30 * 0.70 * 0.43 + 0) = 0.246\n",
      "beta2(2) = sum[a2j * bj(o3) * beta3(j)] = (0.30 * 0.50 * 0.46 + 0.50 * 0.40 * 0.51 + 0.20 * 0.70 * 0.43 + 0) = 0.231\n",
      "beta2(3) = sum[a3j * bj(o3) * beta3(j)] = (0.20 * 0.50 * 0.46 + 0.30 * 0.40 * 0.51 + 0.50 * 0.70 * 0.43 + 0) = 0.258\n",
      "beta1(1) = sum[a1j * bj(o2) * beta2(j)] = (0.50 * 0.50 * 0.25 + 0.20 * 0.60 * 0.23 + 0.30 * 0.30 * 0.26 + 0) = 0.112\n",
      "beta1(2) = sum[a2j * bj(o2) * beta2(j)] = (0.30 * 0.50 * 0.25 + 0.50 * 0.60 * 0.23 + 0.20 * 0.30 * 0.26 + 0) = 0.122\n",
      "beta1(3) = sum[a3j * bj(o2) * beta2(j)] = (0.20 * 0.50 * 0.25 + 0.30 * 0.60 * 0.23 + 0.50 * 0.30 * 0.26 + 0) = 0.105\n",
      "\n",
      "观测序列概率：\n",
      "P(O|lambda) = 0.2 * 0.5 * 0.11246 + 0.4 * 0.4 * 0.12174 + 0.4 * 0.7 * 0.10488 + 0 = [0.0600908]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "class HiddenMarkovBackward:\n",
    "    def __init__(self):\n",
    "        self.betas = None\n",
    "        self.backward_P = None\n",
    "\n",
    "    def backward(self, Q, V, A, B, O, PI):\n",
    "        \"\"\"\n",
    "        后向算法\n",
    "        :param Q: 所有可能的状态集合\n",
    "        :param V: 所有可能的观测集合\n",
    "        :param A: 状态转移概率矩阵\n",
    "        :param B: 观测概率矩阵\n",
    "        :param O: 观测序列\n",
    "        :param PI: 初始状态概率向量\n",
    "        \"\"\"\n",
    "        # 状态序列的大小\n",
    "        N = len(Q)\n",
    "        # 观测序列的大小\n",
    "        M = len(O)\n",
    "        # (1)初始化后向概率beta值，书中第201页公式(10.19)\n",
    "        betas = np.ones((N, M))\n",
    "        self.print_betas_T(N, M)\n",
    "\n",
    "        # (2)对观测序列逆向遍历，M-2即为T-1\n",
    "        print(\"\\n从时刻T-1到1观测序列的后向概率：\")\n",
    "        for t in range(M - 2, -1, -1):\n",
    "            # 得到序列对应的索引\n",
    "            index_of_o = V.index(O[t + 1])\n",
    "            # 遍历状态序列\n",
    "            for i in range(N):\n",
    "                # 书中第201页公式(10.20)\n",
    "                betas[i][t] = np.dot(np.multiply(A[i], [b[index_of_o] for b in B]),\n",
    "                                     [beta[t + 1] for beta in betas])\n",
    "                real_t = t + 1\n",
    "                real_i = i + 1\n",
    "                self.print_betas_t(A, B, N, betas, i, index_of_o, real_i, real_t, t)\n",
    "\n",
    "        # 取出第一个值索引，用于得到o1\n",
    "        index_of_o = V.index(O[0])\n",
    "        self.betas = betas\n",
    "        # 书中第201页公式(10.21)\n",
    "        P = np.dot(np.multiply(PI, [b[index_of_o] for b in B]),\n",
    "                   [beta[0] for beta in betas])\n",
    "        self.backward_P = P\n",
    "        self.print_P(B, N, P, PI, betas, index_of_o)\n",
    "\n",
    "    @staticmethod\n",
    "    def print_P(B, N, P, PI, betas, index_of_o):\n",
    "        print(\"\\n观测序列概率：\")\n",
    "        print(\"P(O|lambda) = \", end=\"\")\n",
    "        for i in range(N):\n",
    "            print(\"%.1f * %.1f * %.5f + \"\n",
    "                  % (PI[0][i], B[i][index_of_o], betas[i][0]), end=\"\")\n",
    "        print(\"0 = {}\".format(P))\n",
    "\n",
    "    @staticmethod\n",
    "    def print_betas_t(A, B, N, betas, i, index_of_o, real_i, real_t, t):\n",
    "        print(\"beta%d(%d) = sum[a%dj * bj(o%d) * beta%d(j)] = (\"\n",
    "              % (real_t, real_i, real_i, real_t + 1, real_t + 1), end='')\n",
    "        for j in range(N):\n",
    "            print(\"%.2f * %.2f * %.2f + \"\n",
    "                  % (A[i][j], B[j][index_of_o], betas[j][t + 1]), end='')\n",
    "        print(\"0) = %.3f\" % betas[i][t])\n",
    "\n",
    "    @staticmethod\n",
    "    def print_betas_T(N, M):\n",
    "        print(\"初始化后向概率：\")\n",
    "        for i in range(N):\n",
    "            print('beta%d(%d) = 1' % (M, i + 1))\n",
    "\n",
    "            \n",
    "Q = [1, 2, 3]\n",
    "V = ['红', '白']\n",
    "A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]\n",
    "B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]\n",
    "O = ['红', '白', '红', '白']\n",
    "PI = [[0.2, 0.4, 0.4]]\n",
    "\n",
    "hmm_backward = HiddenMarkovBackward()\n",
    "hmm_backward.backward(Q, V, A, B, O, PI)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3f880b08",
   "metadata": {},
   "source": [
    "## 习题10.2\n",
    "![image.png](./images/exercise2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5ce2b61b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "前向概率初值：\n",
      "alpha1(1) = pi0 * b0 * b(o1) = 0.100000\n",
      "alpha1(2) = pi1 * b1 * b(o1) = 0.120000\n",
      "alpha1(3) = pi2 * b2 * b(o1) = 0.350000\n",
      "\n",
      "从时刻1到T-1观测序列的前向概率：\n",
      "alpha2(1) = [sum alpha0(i) * ai0] * b0(o1) = 0.078000\n",
      "alpha2(2) = [sum alpha0(i) * ai1] * b1(o1) = 0.084000\n",
      "alpha2(3) = [sum alpha0(i) * ai2] * b2(o1) = 0.082200\n",
      "alpha3(1) = [sum alpha1(i) * ai0] * b0(o2) = 0.040320\n",
      "alpha3(2) = [sum alpha1(i) * ai1] * b1(o2) = 0.026496\n",
      "alpha3(3) = [sum alpha1(i) * ai2] * b2(o2) = 0.068124\n",
      "alpha4(1) = [sum alpha2(i) * ai0] * b0(o3) = 0.020867\n",
      "alpha4(2) = [sum alpha2(i) * ai1] * b1(o3) = 0.012362\n",
      "alpha4(3) = [sum alpha2(i) * ai2] * b2(o3) = 0.043611\n",
      "alpha5(1) = [sum alpha3(i) * ai0] * b0(o4) = 0.011432\n",
      "alpha5(2) = [sum alpha3(i) * ai1] * b1(o4) = 0.010194\n",
      "alpha5(3) = [sum alpha3(i) * ai2] * b2(o4) = 0.011096\n",
      "alpha6(1) = [sum alpha4(i) * ai0] * b0(o5) = 0.005497\n",
      "alpha6(2) = [sum alpha4(i) * ai1] * b1(o5) = 0.003384\n",
      "alpha6(3) = [sum alpha4(i) * ai2] * b2(o5) = 0.009288\n",
      "alpha7(1) = [sum alpha5(i) * ai0] * b0(o6) = 0.002811\n",
      "alpha7(2) = [sum alpha5(i) * ai1] * b1(o6) = 0.002460\n",
      "alpha7(3) = [sum alpha5(i) * ai2] * b2(o6) = 0.002535\n",
      "alpha8(1) = [sum alpha6(i) * ai0] * b0(o7) = 0.001325\n",
      "alpha8(2) = [sum alpha6(i) * ai1] * b1(o7) = 0.001211\n",
      "alpha8(3) = [sum alpha6(i) * ai2] * b2(o7) = 0.000941\n",
      "\n",
      "初始化后向概率：\n",
      "beta8(1) = 1\n",
      "beta8(2) = 1\n",
      "beta8(3) = 1\n",
      "\n",
      "从时刻T-1到1观测序列的后向概率：\n",
      "beta7(1) = sum[a1j * bj(o8) * beta8(j)] = (0.50 * 0.50 * 1.00 + 0.10 * 0.60 * 1.00 + 0.40 * 0.30 * 1.00 + 0) = 0.430\n",
      "beta7(2) = sum[a2j * bj(o8) * beta8(j)] = (0.30 * 0.50 * 1.00 + 0.50 * 0.60 * 1.00 + 0.20 * 0.30 * 1.00 + 0) = 0.510\n",
      "beta7(3) = sum[a3j * bj(o8) * beta8(j)] = (0.20 * 0.50 * 1.00 + 0.20 * 0.60 * 1.00 + 0.60 * 0.30 * 1.00 + 0) = 0.400\n",
      "beta6(1) = sum[a1j * bj(o7) * beta7(j)] = (0.50 * 0.50 * 0.43 + 0.10 * 0.60 * 0.51 + 0.40 * 0.30 * 0.40 + 0) = 0.186\n",
      "beta6(2) = sum[a2j * bj(o7) * beta7(j)] = (0.30 * 0.50 * 0.43 + 0.50 * 0.60 * 0.51 + 0.20 * 0.30 * 0.40 + 0) = 0.241\n",
      "beta6(3) = sum[a3j * bj(o7) * beta7(j)] = (0.20 * 0.50 * 0.43 + 0.20 * 0.60 * 0.51 + 0.60 * 0.30 * 0.40 + 0) = 0.176\n",
      "beta5(1) = sum[a1j * bj(o6) * beta6(j)] = (0.50 * 0.50 * 0.19 + 0.10 * 0.40 * 0.24 + 0.40 * 0.70 * 0.18 + 0) = 0.106\n",
      "beta5(2) = sum[a2j * bj(o6) * beta6(j)] = (0.30 * 0.50 * 0.19 + 0.50 * 0.40 * 0.24 + 0.20 * 0.70 * 0.18 + 0) = 0.101\n",
      "beta5(3) = sum[a3j * bj(o6) * beta6(j)] = (0.20 * 0.50 * 0.19 + 0.20 * 0.40 * 0.24 + 0.60 * 0.70 * 0.18 + 0) = 0.112\n",
      "beta4(1) = sum[a1j * bj(o5) * beta5(j)] = (0.50 * 0.50 * 0.11 + 0.10 * 0.60 * 0.10 + 0.40 * 0.30 * 0.11 + 0) = 0.046\n",
      "beta4(2) = sum[a2j * bj(o5) * beta5(j)] = (0.30 * 0.50 * 0.11 + 0.50 * 0.60 * 0.10 + 0.20 * 0.30 * 0.11 + 0) = 0.053\n",
      "beta4(3) = sum[a3j * bj(o5) * beta5(j)] = (0.20 * 0.50 * 0.11 + 0.20 * 0.60 * 0.10 + 0.60 * 0.30 * 0.11 + 0) = 0.043\n",
      "beta3(1) = sum[a1j * bj(o4) * beta4(j)] = (0.50 * 0.50 * 0.05 + 0.10 * 0.40 * 0.05 + 0.40 * 0.70 * 0.04 + 0) = 0.026\n",
      "beta3(2) = sum[a2j * bj(o4) * beta4(j)] = (0.30 * 0.50 * 0.05 + 0.50 * 0.40 * 0.05 + 0.20 * 0.70 * 0.04 + 0) = 0.023\n",
      "beta3(3) = sum[a3j * bj(o4) * beta4(j)] = (0.20 * 0.50 * 0.05 + 0.20 * 0.40 * 0.05 + 0.60 * 0.70 * 0.04 + 0) = 0.027\n",
      "beta2(1) = sum[a1j * bj(o3) * beta3(j)] = (0.50 * 0.50 * 0.03 + 0.10 * 0.40 * 0.02 + 0.40 * 0.70 * 0.03 + 0) = 0.015\n",
      "beta2(2) = sum[a2j * bj(o3) * beta3(j)] = (0.30 * 0.50 * 0.03 + 0.50 * 0.40 * 0.02 + 0.20 * 0.70 * 0.03 + 0) = 0.012\n",
      "beta2(3) = sum[a3j * bj(o3) * beta3(j)] = (0.20 * 0.50 * 0.03 + 0.20 * 0.40 * 0.02 + 0.60 * 0.70 * 0.03 + 0) = 0.016\n",
      "beta1(1) = sum[a1j * bj(o2) * beta2(j)] = (0.50 * 0.50 * 0.01 + 0.10 * 0.60 * 0.01 + 0.40 * 0.30 * 0.02 + 0) = 0.006\n",
      "beta1(2) = sum[a2j * bj(o2) * beta2(j)] = (0.30 * 0.50 * 0.01 + 0.50 * 0.60 * 0.01 + 0.20 * 0.30 * 0.02 + 0) = 0.007\n",
      "beta1(3) = sum[a3j * bj(o2) * beta2(j)] = (0.20 * 0.50 * 0.01 + 0.20 * 0.60 * 0.01 + 0.60 * 0.30 * 0.02 + 0) = 0.006\n",
      "\n",
      "观测序列概率：\n",
      "P(O|lambda) = 0.2 * 0.5 * 0.00633 + 0.3 * 0.4 * 0.00685 + 0.5 * 0.7 * 0.00578 + 0 = [0.00347671]\n",
      "\n",
      "计算P(i4=q3|O,lambda)：\n",
      "P(i4=q3|O,lambda) = alpha4(3) * beta4(3) / P(O|lambda) = 0.536952\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "class HiddenMarkovForwardBackward(HiddenMarkovBackward):\n",
    "    def __init__(self, verbose=False):\n",
    "        super(HiddenMarkovBackward, self).__init__()\n",
    "        self.alphas = None\n",
    "        self.forward_P = None\n",
    "        self.verbose = verbose\n",
    "\n",
    "    def forward(self, Q, V, A, B, O, PI):\n",
    "        \"\"\"\n",
    "        前向算法\n",
    "        :param Q: 所有可能的状态集合\n",
    "        :param V: 所有可能的观测集合\n",
    "        :param A: 状态转移概率矩阵\n",
    "        :param B: 观测概率矩阵\n",
    "        :param O: 观测序列\n",
    "        :param PI: 初始状态概率向量\n",
    "        \"\"\"\n",
    "        # 状态序列的大小\n",
    "        N = len(Q)\n",
    "        # 观测序列的大小\n",
    "        M = len(O)\n",
    "        # 初始化前向概率alpha值\n",
    "        alphas = np.zeros((N, M))\n",
    "        # 时刻数=观测序列数\n",
    "        T = M\n",
    "        # (2)对观测序列遍历，遍历每一个时刻，计算前向概率alpha值\n",
    "\n",
    "        for t in range(T):\n",
    "            if self.verbose:\n",
    "                if t == 0:\n",
    "                    print(\"前向概率初值：\")\n",
    "                elif t == 1:\n",
    "                    print(\"\\n从时刻1到T-1观测序列的前向概率：\")\n",
    "            # 得到序列对应的索引\n",
    "            index_of_o = V.index(O[t])\n",
    "            # 遍历状态序列\n",
    "            for i in range(N):\n",
    "                if t == 0:\n",
    "                    # (1)初始化alpha初值，书中第198页公式(10.15)\n",
    "                    alphas[i][t] = PI[t][i] * B[i][index_of_o]\n",
    "                    if self.verbose:\n",
    "                        self.print_alpha_t1(alphas, i, t)\n",
    "                else:\n",
    "                    # (2)递推，书中第198页公式(10.16)\n",
    "                    alphas[i][t] = np.dot([alpha[t - 1] for alpha in alphas],\n",
    "                                          [a[i] for a in A]) * B[i][index_of_o]\n",
    "                    if self.verbose:\n",
    "                        self.print_alpha_t(alphas, i, t)\n",
    "        # (3)终止，书中第198页公式(10.17)\n",
    "        self.forward_P = np.sum([alpha[M - 1] for alpha in alphas])\n",
    "        self.alphas = alphas\n",
    "\n",
    "    @staticmethod\n",
    "    def print_alpha_t(alphas, i, t):\n",
    "        print(\"alpha%d(%d) = [sum alpha%d(i) * ai%d] * b%d(o%d) = %f\"\n",
    "              % (t + 1, i + 1, t - 1, i, i, t, alphas[i][t]))\n",
    "\n",
    "    @staticmethod\n",
    "    def print_alpha_t1(alphas, i, t):\n",
    "        print('alpha1(%d) = pi%d * b%d * b(o1) = %f'\n",
    "              % (i + 1, i, i, alphas[i][t]))\n",
    "\n",
    "    def calc_t_qi_prob(self, t, qi):\n",
    "        result = (self.alphas[qi - 1][t - 1] * self.betas[qi - 1][t - 1]) / self.backward_P[0]\n",
    "        if self.verbose:\n",
    "            print(\"计算P(i%d=q%d|O,lambda)：\" % (t, qi))\n",
    "            print(\"P(i%d=q%d|O,lambda) = alpha%d(%d) * beta%d(%d) / P(O|lambda) = %f\"\n",
    "                  % (t, qi, t, qi, t, qi, result))\n",
    "\n",
    "        return result\n",
    "\n",
    "    \n",
    "Q = [1, 2, 3]\n",
    "V = ['红', '白']\n",
    "A = [[0.5, 0.1, 0.4], [0.3, 0.5, 0.2], [0.2, 0.2, 0.6]]\n",
    "B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]\n",
    "O = ['红', '白', '红', '红', '白', '红', '白', '白']\n",
    "PI = [[0.2, 0.3, 0.5]]\n",
    "\n",
    "hmm_forward_backward = HiddenMarkovForwardBackward(verbose=True)\n",
    "hmm_forward_backward.forward(Q, V, A, B, O, PI)\n",
    "print()\n",
    "hmm_forward_backward.backward(Q, V, A, B, O, PI)\n",
    "print()\n",
    "hmm_forward_backward.calc_t_qi_prob(t=4, qi=3)\n",
    "print()\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "70e3399b",
   "metadata": {},
   "source": [
    "## 习题10.3\n",
    "![image.png](./images/exercise3.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c46b3542",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "初始化Psi1和delta1：\n",
      "delta1(1) = pi1 * b1(o1) = 0.20 * 0.50 = 0.10\n",
      "Psi1(1) = 0\n",
      "delta1(2) = pi2 * b2(o1) = 0.40 * 0.40 = 0.16\n",
      "Psi1(2) = 0\n",
      "delta1(3) = pi3 * b3(o1) = 0.40 * 0.70 = 0.28\n",
      "Psi1(3) = 0\n",
      "\n",
      "从时刻2到T的所有单个路径中概率最大值delta和概率最大的路径的第t-1个结点Psi：\n",
      "delta2(1) = max[delta1(j) * aj1] * b1(o2) = 0.06 * 0.50 = 0.02800\n",
      "Psi2(1) = argmax[delta1(j) * aj1] = 2\n",
      "delta2(2) = max[delta1(j) * aj2] * b2(o2) = 0.08 * 0.60 = 0.05040\n",
      "Psi2(2) = argmax[delta1(j) * aj2] = 2\n",
      "delta2(3) = max[delta1(j) * aj3] * b3(o2) = 0.14 * 0.30 = 0.04200\n",
      "Psi2(3) = argmax[delta1(j) * aj3] = 2\n",
      "delta3(1) = max[delta2(j) * aj1] * b1(o3) = 0.02 * 0.50 = 0.00756\n",
      "Psi3(1) = argmax[delta2(j) * aj1] = 1\n",
      "delta3(2) = max[delta2(j) * aj2] * b2(o3) = 0.03 * 0.40 = 0.01008\n",
      "Psi3(2) = argmax[delta2(j) * aj2] = 1\n",
      "delta3(3) = max[delta2(j) * aj3] * b3(o3) = 0.02 * 0.70 = 0.01470\n",
      "Psi3(3) = argmax[delta2(j) * aj3] = 2\n",
      "delta4(1) = max[delta3(j) * aj1] * b1(o4) = 0.00 * 0.50 = 0.00189\n",
      "Psi4(1) = argmax[delta3(j) * aj1] = 0\n",
      "delta4(2) = max[delta3(j) * aj2] * b2(o4) = 0.01 * 0.60 = 0.00302\n",
      "Psi4(2) = argmax[delta3(j) * aj2] = 1\n",
      "delta4(3) = max[delta3(j) * aj3] * b3(o4) = 0.01 * 0.30 = 0.00220\n",
      "Psi4(3) = argmax[delta3(j) * aj3] = 2\n",
      "\n",
      "所有路径的终结点最大的概率值：\n",
      "P = 0.003024\n",
      "\n",
      "最优路径的终结点：\n",
      "i4 = argmax[deltaT(i)] = 2\n",
      "\n",
      "最优路径的其他结点：\n",
      "i3 = Psi4(i4) = 2\n",
      "i2 = Psi3(i3) = 2\n",
      "i1 = Psi2(i2) = 3\n",
      "\n",
      "最优路径是： 3->2->2->2\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "class HiddenMarkovViterbi:\n",
    "    def __init__(self, verbose=False):\n",
    "        self.verbose = verbose\n",
    "\n",
    "    def viterbi(self, Q, V, A, B, O, PI):\n",
    "        \"\"\"\n",
    "        维特比算法\n",
    "        :param Q: 所有可能的状态集合\n",
    "        :param V: 所有可能的观测集合\n",
    "        :param A: 状态转移概率矩阵\n",
    "        :param B: 观测概率矩阵\n",
    "        :param O: 观测序列\n",
    "        :param PI: 初始状态概率向量\n",
    "        \"\"\"\n",
    "        # 状态序列的大小\n",
    "        N = len(Q)\n",
    "        # 观测序列的大小\n",
    "        M = len(O)\n",
    "        # 初始化deltas\n",
    "        deltas = np.zeros((N, M))\n",
    "        # 初始化psis\n",
    "        psis = np.zeros((N, M))\n",
    "\n",
    "        # 初始化最优路径矩阵，该矩阵维度与观测序列维度相同\n",
    "        I = np.zeros((1, M))\n",
    "        # (2)递推，遍历观测序列\n",
    "        for t in range(M):\n",
    "            if self.verbose:\n",
    "                if t == 0:\n",
    "                    print(\"初始化Psi1和delta1：\")\n",
    "                elif t == 1:\n",
    "                    print(\"\\n从时刻2到T的所有单个路径中概率\"\n",
    "                          \"最大值delta和概率最大的路径的第t-1个结点Psi：\")\n",
    "\n",
    "            # (2)递推从t=2开始\n",
    "            real_t = t + 1\n",
    "            # 得到序列对应的索引\n",
    "            index_of_o = V.index(O[t])\n",
    "            for i in range(N):\n",
    "                real_i = i + 1\n",
    "                if t == 0:\n",
    "                    # (1)初始化\n",
    "                    deltas[i][t] = PI[0][i] * B[i][index_of_o]\n",
    "                    psis[i][t] = 0\n",
    "\n",
    "                    self.print_delta_t1(\n",
    "                        B, PI, deltas, i, index_of_o, real_i, t)\n",
    "                    self.print_psi_t1(real_i)\n",
    "                else:\n",
    "                    # (2)递推，对t=2,3,...,T\n",
    "                    deltas[i][t] = np.max(np.multiply([delta[t - 1] for delta in deltas],\n",
    "                                                      [a[i] for a in A])) * B[i][index_of_o]\n",
    "                    self.print_delta_t(\n",
    "                        A, B, deltas, i, index_of_o, real_i, real_t, t)\n",
    "\n",
    "                    psis[i][t] = np.argmax(np.multiply([delta[t - 1] for delta in deltas],\n",
    "                                                       [a[i] for a in A]))\n",
    "                    self.print_psi_t(i, psis, real_i, real_t, t)\n",
    "\n",
    "        last_deltas = [delta[M - 1] for delta in deltas]\n",
    "        # (3)终止，得到所有路径的终结点最大的概率值\n",
    "        P = np.max(last_deltas)\n",
    "        # (3)得到最优路径的终结点\n",
    "        I[0][M - 1] = np.argmax(last_deltas)\n",
    "        if self.verbose:\n",
    "            print(\"\\n所有路径的终结点最大的概率值：\")\n",
    "            print(\"P = %f\" % P)\n",
    "        if self.verbose:\n",
    "            print(\"\\n最优路径的终结点：\")\n",
    "            print(\"i%d = argmax[deltaT(i)] = %d\" % (M, I[0][M - 1] + 1))\n",
    "            print(\"\\n最优路径的其他结点：\")\n",
    "\n",
    "        # (4)递归由后向前得到其他结点\n",
    "        for t in range(M - 2, -1, -1):\n",
    "            I[0][t] = psis[int(I[0][t + 1])][t + 1]\n",
    "            if self.verbose:\n",
    "                print(\"i%d = Psi%d(i%d) = %d\" %\n",
    "                      (t + 1, t + 2, t + 2, I[0][t] + 1))\n",
    "\n",
    "        # 输出最优路径\n",
    "        print(\"\\n最优路径是：\", \"->\".join([str(int(i + 1)) for i in I[0]]))\n",
    "\n",
    "    def print_psi_t(self, i, psis, real_i, real_t, t):\n",
    "        if self.verbose:\n",
    "            print(\"Psi%d(%d) = argmax[delta%d(j) * aj%d] = %d\"\n",
    "                  % (real_t, real_i, real_t - 1, real_i, psis[i][t]))\n",
    "\n",
    "    def print_delta_t(self, A, B, deltas, i, index_of_o, real_i, real_t, t):\n",
    "        if self.verbose:\n",
    "            print(\"delta%d(%d) = max[delta%d(j) * aj%d] * b%d(o%d) = %.2f * %.2f = %.5f\"\n",
    "                  % (real_t, real_i, real_t - 1, real_i, real_i, real_t,\n",
    "                     np.max(np.multiply([delta[t - 1] for delta in deltas],\n",
    "                                        [a[i] for a in A])),\n",
    "                     B[i][index_of_o], deltas[i][t]))\n",
    "\n",
    "    def print_psi_t1(self, real_i):\n",
    "        if self.verbose:\n",
    "            print(\"Psi1(%d) = 0\" % real_i)\n",
    "\n",
    "    def print_delta_t1(self, B, PI, deltas, i, index_of_o, real_i, t):\n",
    "        if self.verbose:\n",
    "            print(\"delta1(%d) = pi%d * b%d(o1) = %.2f * %.2f = %.2f\"\n",
    "                  % (real_i, real_i, real_i, PI[0][i], B[i][index_of_o], deltas[i][t]))\n",
    "\n",
    "            \n",
    "            \n",
    "Q = [1, 2, 3]\n",
    "V = ['红', '白']\n",
    "A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]\n",
    "B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]\n",
    "O = ['红', '白', '红', '白']\n",
    "PI = [[0.2, 0.4, 0.4]]\n",
    "\n",
    "HMM = HiddenMarkovViterbi(verbose=True)\n",
    "HMM.viterbi(Q, V, A, B, O, PI)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23e8c082",
   "metadata": {},
   "source": [
    "# 自编程实现隐马尔可夫模型的前向、后向、维特比算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "39be6edd",
   "metadata": {},
   "outputs": [],
   "source": [
    "class HiddenMarkov:\n",
    "    def __init__(self):\n",
    "        self.alphas = None\n",
    "        self.forward_P = None\n",
    "        self.betas = None\n",
    "        self.backward_P = None\n",
    "\n",
    "    # 前向算法\n",
    "    def forward(self, Q, V, A, B, O, PI):\n",
    "        # 状态序列的大小\n",
    "        N = len(Q)\n",
    "        # 观测序列的大小\n",
    "        M = len(O)\n",
    "        # 初始化前向概率alpha值\n",
    "        alphas = np.zeros((N, M))\n",
    "        # 时刻数=观测序列数\n",
    "        T = M\n",
    "        # 遍历每一个时刻，计算前向概率alpha值\n",
    "        for t in range(T):\n",
    "            # 得到序列对应的索引\n",
    "            indexOfO = V.index(O[t])\n",
    "            # 遍历状态序列\n",
    "            for i in range(N):\n",
    "                # 初始化alpha初值\n",
    "                if t == 0:\n",
    "                    # P176 公式(10.15)\n",
    "                    alphas[i][t] = PI[t][i] * B[i][indexOfO]\n",
    "                    print('alpha1(%d) = p%db%db(o1) = %f' %\n",
    "                          (i + 1, i, i, alphas[i][t]))\n",
    "                else:\n",
    "                    # P176 公式(10.16)\n",
    "                    alphas[i][t] = np.dot([alpha[t - 1] for alpha in alphas],\n",
    "                                          [a[i] for a in A]) * B[i][indexOfO]\n",
    "                    print('alpha%d(%d) = [sigma alpha%d(i)ai%d]b%d(o%d) = %f' %\n",
    "                          (t + 1, i + 1, t - 1, i, i, t, alphas[i][t]))\n",
    "        # P176 公式(10.17)\n",
    "        self.forward_P = np.sum([alpha[M - 1] for alpha in alphas])\n",
    "        self.alphas = alphas\n",
    "\n",
    "    # 后向算法\n",
    "    def backward(self, Q, V, A, B, O, PI):\n",
    "        # 状态序列的大小\n",
    "        N = len(Q)\n",
    "        # 观测序列的大小\n",
    "        M = len(O)\n",
    "        # 初始化后向概率beta值，P178 公式(10.19)\n",
    "        betas = np.ones((N, M))\n",
    "        #\n",
    "        for i in range(N):\n",
    "            print('beta%d(%d) = 1' % (M, i + 1))\n",
    "        # 对观测序列逆向遍历\n",
    "        for t in range(M - 2, -1, -1):\n",
    "            # 得到序列对应的索引\n",
    "            indexOfO = V.index(O[t + 1])\n",
    "            # 遍历状态序列\n",
    "            for i in range(N):\n",
    "                # P178 公式(10.20)\n",
    "                betas[i][t] = np.dot(\n",
    "                    np.multiply(A[i], [b[indexOfO] for b in B]),\n",
    "                    [beta[t + 1] for beta in betas])\n",
    "                realT = t + 1\n",
    "                realI = i + 1\n",
    "                print('beta%d(%d) = sigma[a%djbj(o%d)beta%d(j)] = (' %\n",
    "                      (realT, realI, realI, realT + 1, realT + 1),\n",
    "                      end='')\n",
    "                for j in range(N):\n",
    "                    print(\"%.2f * %.2f * %.2f + \" %\n",
    "                          (A[i][j], B[j][indexOfO], betas[j][t + 1]),\n",
    "                          end='')\n",
    "                print(\"0) = %.3f\" % betas[i][t])\n",
    "        # 取出第一个值\n",
    "        indexOfO = V.index(O[0])\n",
    "        self.betas = betas\n",
    "        # P178 公式(10.21)\n",
    "        P = np.dot(np.multiply(PI, [b[indexOfO] for b in B]),\n",
    "                   [beta[0] for beta in betas])\n",
    "        self.backward_P = P\n",
    "        print(\"P(O|lambda) = \", end=\"\")\n",
    "        for i in range(N):\n",
    "            print(\"%.1f * %.1f * %.5f + \" %\n",
    "                  (PI[0][i], B[i][indexOfO], betas[i][0]),\n",
    "                  end=\"\")\n",
    "        print(\"0 = %f\" % P)\n",
    "\n",
    "    # 维特比算法\n",
    "    def viterbi(self, Q, V, A, B, O, PI):\n",
    "        # 状态序列的大小\n",
    "        N = len(Q)\n",
    "        # 观测序列的大小\n",
    "        M = len(O)\n",
    "        # 初始化daltas\n",
    "        deltas = np.zeros((N, M))\n",
    "        # 初始化psis\n",
    "        psis = np.zeros((N, M))\n",
    "        # 初始化最优路径矩阵，该矩阵维度与观测序列维度相同\n",
    "        I = np.zeros((1, M))\n",
    "        # 遍历观测序列\n",
    "        for t in range(M):\n",
    "            # 递推从t=2开始\n",
    "            realT = t + 1\n",
    "            # 得到序列对应的索引\n",
    "            indexOfO = V.index(O[t])\n",
    "            for i in range(N):\n",
    "                realI = i + 1\n",
    "                if t == 0:\n",
    "                    # P185 算法10.5 步骤(1)\n",
    "                    deltas[i][t] = PI[0][i] * B[i][indexOfO]\n",
    "                    psis[i][t] = 0\n",
    "                    print('delta1(%d) = pi%d * b%d(o1) = %.2f * %.2f = %.2f' %\n",
    "                          (realI, realI, realI, PI[0][i], B[i][indexOfO],\n",
    "                           deltas[i][t]))\n",
    "                    print('psis1(%d) = 0' % (realI))\n",
    "                else:\n",
    "                    # # P185 算法10.5 步骤(2)\n",
    "                    deltas[i][t] = np.max(\n",
    "                        np.multiply([delta[t - 1] for delta in deltas],\n",
    "                                    [a[i] for a in A])) * B[i][indexOfO]\n",
    "                    print(\n",
    "                        'delta%d(%d) = max[delta%d(j)aj%d]b%d(o%d) = %.2f * %.2f = %.5f'\n",
    "                        % (realT, realI, realT - 1, realI, realI, realT,\n",
    "                           np.max(\n",
    "                               np.multiply([delta[t - 1] for delta in deltas],\n",
    "                                           [a[i] for a in A])), B[i][indexOfO],\n",
    "                           deltas[i][t]))\n",
    "                    psis[i][t] = np.argmax(\n",
    "                        np.multiply([delta[t - 1] for delta in deltas],\n",
    "                                    [a[i] for a in A]))\n",
    "                    print('psis%d(%d) = argmax[delta%d(j)aj%d] = %d' %\n",
    "                          (realT, realI, realT - 1, realI, psis[i][t]))\n",
    "        #print(deltas)\n",
    "        #print(psis)\n",
    "        # 得到最优路径的终结点\n",
    "        I[0][M - 1] = np.argmax([delta[M - 1] for delta in deltas])\n",
    "        print('i%d = argmax[deltaT(i)] = %d' % (M, I[0][M - 1] + 1))\n",
    "        # 递归由后向前得到其他结点\n",
    "        for t in range(M - 2, -1, -1):\n",
    "            I[0][t] = psis[int(I[0][t + 1])][t + 1]\n",
    "            print('i%d = psis%d(i%d) = %d' %\n",
    "                  (t + 1, t + 2, t + 2, I[0][t] + 1))\n",
    "        # 输出最优路径\n",
    "        print('最优路径是：', \"->\".join([str(int(i + 1)) for i in I[0]]))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
